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Abstract. Determining the equilibrium configuration of an elastic Mobius band is a 
challenging problem. In recent years numerical results have been obtained by other 
investigators, employing first the Kirchhoff theory of rods and later the developable, ruled- 
surface model of Wunderlich. In pailicular, the strategy employed previously for the latter does 
not deliver an unconstrained equilibrium configuration for the complete strip. Here we present 
our own systematic approach to the same problem for each of these models, with the ultimate 
goal of assessing the stability of flip-symmetric configurations. The presence of pointwise 
constraints considerably complicates the latter step. We obtain the first stability results for 
the problem, numerically demonstrating that such equilibria render the total potential energy 
a local minimum. Along the way we introduce a novel regularization for the for the singular 
Wunderlich model that delivers unconstrained equilibria for the complete strip, which can then 
be tested for stability. 


1. Introduction 

Sadowsky |[T]|2l and later Wunderlich jS) were the hrst to propose models for determining the 
equilibrium conhgurations of elastic Mobius bands, idealizing them as a developable surfaces. 
In particular, linear isotropic plate theory is employed in El where an integration across the 
width yields an energy density per unit length, reminiscent of rod theory. Obtaining numerical 
solutions for the latter is challenging and was only taken up recently, cf ||4l, El- The 
interpretation of the model derived in El in light of classical Kirchhoff rod theory was made 
precise in the recent work |0. In an earlier work 17], Kirchhoff rod theory was employed to 
obtain certain smooth shapes of Mobius strips, in which one cross-sectional moment of inertia 
of the rod is much larger than the other. 

In both Q and m, El the equilibrium equations are solved numerically on the half¬ 
domain with appropriate boundary conditions; assuming flip symmetry, the full solution is 
generated by rotation through n radians about the symmetry axis. That is, the entire closed- 
loop configuration has a single rotational symmetry (by n radians) about some fixed axis. 
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The same approach was employed to obtain conhgurations of twisted isotropic rods in IS]. 
There the two ends of a straight isotropic rod are twisted through any relative angle and then 
seamlessly joined, as compared to the Mobius band, where the relative angle is necessarily 
7Z radians. In fS] it is shown that all equilibrium conhgurations possess hip symmetry, i.e., 
the above solution procedure can be used without loss of generality. That argument relies 
crucially on cross-sectional isotropy. In particular, we know of no such result here for strips, 
i.e., other non-symmetrical solutions could exist. 

Given that, we address a more modest but nonetheless important question here, viz., we 
assess the local stability of the hip-symmetric solution. Stable or not, this does not rule out the 
possibility of non-symmetrical solutions. We consider the two distinct models employed in 
|7l and |@), viz., the Kirchhoff rod model and the Wunderlich model, respectively. The former 
serves as a “warm-up” for the latter. Also the Kirchhoff model is a reasonable one for bands 
made of compliant materials like rubber. In order to test stability, we must hrst obtain reliable 
symmetric equilibria for the full closed loop, which is a challenging task. This is particularly 
true for the developable-surface model, due to the inherent singularity associated with the 
ruled-surface parametrization employed in 0. For various reasons, discussed below, we do 
not employ the formulations of |2l and 0], |5]. Accordingly, the paper is taken up presenting 
our systematic formulations for both the numerical computation of symmetric equilibria and 
the assessment of their stability. While the latter is certainly new, the former constitutes the 
hrst systematic approach to computing unconstrained equilibria of complete strips for the 
developable-surface model of Wunderlich. 

The outline of this work is as follows. In Section 2 we summarize the well-known held 
equations for hyperelastic, inextensible, unshearable Cosserat rods, ultimately adopting the 
constitutive assumption normally attributed to Kirchhoff. In Section 3 we summarize our 
formulation, as hrst presented in |9l, and compute solutions for the half rod with appropriate 
boundary conditions engendering hip symmetry. We avoid the use of Euler angles and their 
associated singularities as in |2] and |5l; the kinematical description of the hnite rotation held 
here is singularity-free via quaternions. As shown in |9l, the exploitation of a “conserved” 
quantity delivers a complete formulation within the context of a linear space. As in jT] we use 
the ratio of the cross-sectional area moments of inertia as a continuation parameter - starting 
from the well-known hat, circular equilibrium conhguration. In anticipation of our stability 
results, we then extend all solution helds - kinematic and kinetic - to the entire closed-loop 
conhguration. In Section 4 we briehy present our results for hip-symmetric equilibria. 

In a conservative problem such as the one at hand, it’s enough to check the positivity 
of the (reduced) stiffness matrix at an equilibrium to deduce that the total potential energy is 
a local minimum there, i.e., the conhguration is locally stable. Unfortunately in the case of 
two-point boundary value problems, such information is not a direct by-product of the code 
AUTO. However, the real difficulty here stems from pointwise constraints like inextensibility 
and unshearability, present in the problem at hand. In Section 5 we employ the methodology 
of ifTOl to overcome this. We hrst identify the discrete, numerical solution for the closed 
loop with a hnite-element mesh, and then consider its linearization about the equilibrium 
conhguration. This yields a stiffness matrix in the presence of constraints. A QR-factorization 
of the constraint matrix enables the determination of the symmetric projected stiffness matrix, 
dehned on the orthogonal complement of the subspace spanned by the constraints. We then 
compute the smallest eigenvalues of the projected stiffness matrix. In this way we numerically 
verify the stability of all closed-loop solutions found. 

In Section 6 we take up the Wunderlich model, with the same goals in mind as above. As 
hrst noted in |4i, but more clearly illuminated in |6|, the resulting held equations are those of 
a Cosserat rod in the presence of an additional “state variable” of a purely geometric nature. In 
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(a) Half-Strip Solution (b) Generated Full Strip 


Figure 1. The half-strip solution and the constructed full strip for the developable rod. 
Indicated are the moments niappjied ImappUed in the half strip and full strip respectively 
caused by specifying a non-zero curvature on the half strip’s end. These are external applied 
moments concentrated at the end point of the half strip. 


particular, the governing equation associated with the latter possesses a singularity wherever 
the curvature of the centerline curve vanishes im. As observed in Q, such a condition 
necessarily occurs at one end of the half-rod on the symmetry axis. This renders the numerical 
determination of complete flip-symmetric configurations and their stability assessment much 
more difficult. In Q, a small external curvamre is imposed at one end of the half-band in 
order to overcome the singularity that is otherwise present at that location. This is equivalent 
to the presence of a small externally applied moment at that end. Consequently, this method 
does not yield an unconstrained equilibrium configuration: When the half configuration is 
rotated about the symmetry axis through an angle of n, the small applied moment is doubled 
in magnitude and acts externally on the complete band, as shown in Figure[T] In the structural 
mechanics literature this is often referred to as anti-symmetry, cf. m. In any case, if the 
boundary supports for this half-rod are removed, the full band does not satisfy global moment 
balance. Here we take a different approach to obtain unconstrained equilibria for the full 
strip. We use the same rod formulation described above but now add an internal “elliptic 
regularization” associated with the state variable, characterized by a small parameter, to the 
potential energy derived in iB). The resulting Euler-Lagrange equation for the state variable 
is now singularly perturbed but not singular, and the extension of the solution to a complete 
configuration does not give rise to an unbalanced external moment. With this in hand, we 
carry out the same strategy described above for the Kirchhoff rod, but accompanied by taking 
the regularizing parameter as small as possible in the continuation scheme. 

Another point of departure from the strategy used in the Kirchhoff model is that a known 
configuration to initiate continuation here is not at all obvious. In Section 7 we start from 
a circularly bent, untwisted half strip of fixed width; this requires the application of an end 
moment at the hinged end. We then execute a two-parameter continuation - first relaxing the 
end moment and then aligning the hinge at the appropriate ;r/2 orientation. Hereafter, the 
width and small regularizing parameter are employed as continuation parameters to obtain 
half-band configurations. In Section 8 we present our results for flip-symmetric closed-loop 
configurations. In Section 9 we take up the assessment of their stability. We first extend the 
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development in ifTOl to the more complicated problem at hand. Then after a careful extension 
of all computed fields on the half strip to the full closed loop, the computational procedure for 
obtaining the projected stiffness matrix is the same as above. Again, we find that all computed 
flip-symmetric configurations correspond to local energy minima. 

2. Elastic Rod Formulation 

Let { 61 , 62 , 63 } denote a fixed, right-handed, orthonormal basis for E^, the translate space 
for 3-dimensional Euclidean point space. We start by defining the special Cosserat rod 
with centerline coordinate s G [0,L] in a straight, stress-free reference configuration. The 
position of the rod is defined by the vector-valued function r : [ 0 ,L] —>■ with the reference 
configuration’s centerline given by Tq (i) = ^ 63 . The cross-sections of the rod in the reference 
configuration are parallel to the plane span { 61 , 62 }. Let R (s) denote the rotation of the cross- 
sectional plane parallel to span { 61 , 62 } at i in the undeformed rod. 

We define an orthonormal basis field jd 1 , d 2 , d 3 } via 

d,-(s)=R(s) 6 ,. (1) 

The configuration of the rod is uniquely determined by the functions r(i) and R(i). The 
director basis field, {di ,d 2 ,d 3 } is attached to the centroid of the rod’s cross section and defines 
the orientation of the rod’s cross sections. Differentiation of Q yields 

d; = tfxd,-, ( 2 ) 

where a x b denotes the usual right-handed cross product, / denotes a derivative with respect 
to the centerline coordinate, s, and K is the axial vector field of the skew symmetric tensor 
field R'R^, denoted 

K-= axial (R'R^) . (3) 

We consider here, unshearable and inextensible rods, viz., we impose the constraint 

r' = d 3 . (4) 

We also write 

K=Kidi, (5) 

where here and throughout repeated Latin subscripts imply summation from 1 to 3, while 
repeated Greek subscripts sum from 1 to 2. The scalar fields tc,, i = 1,2,3, are the strains of 
the theory: Ki , K 2 are components of the curvature or bending strains, while Kj is the twist or 
torsional strain. 

The vector fields n (s) and m (s), denote the internal contact force and contact couple, 
respectively, acting on the deformed cross section at “s”. We write 

n = n,d,-, ( 6 ) 

m = m; d;, (7) 

where the component fields, n, and ot;, i = 1,2,3, are the internal forces and moments 
respectively: ni,n 2 correspond to shear forces; n^, axial force; 1111,1112 correspond to bending 


Computation of Unconstrained Elastic Equilibria of Complete Mobius Bands and their StabilityS 


moments, and m^ to torque or twisting moment. In the absence of body forces and body 
couples, the local form of force and moment balance are given by 

n' = 0, (8) 

m' + da X n = 0, (9) 

respectively where we have used Qua. 

Since the tangent vector of the centerline, r' is constrained to be da, the rod contact force, 
n, is not constitutively determined. In other words, the contact force serves as a Lagrange 
multiplier enforcing the unshearable-inextensible constraint 0 - 

We define an objective, hyperelastic, inextensible and unshearable rod as one 
characterized by the existence of a non-negative function W : —> [ 0 , oo), called the stored 
energy density, such that 

dW 

mi=— , i= 1,2,3. (10) 

oKi 

For notational convenience, we denote the following triples of real number via 


k := (K-i,K- 2 ,K-a) n := (ni,n2,n3) m := {m\,m 2 ,mf) 


( 11 ) 


( 12 ) 


Writing W {V.)=W (tfi, K 2 , Kf), then (10 1 takes the compact form 

dW 
dlT' 

With the aid of Q-Q, 0-0 and ( [T^ , the balance equations (|^ and Q take the form 

n' -f k X n = 0, (13) 

m'+ k X m-fd X n = 0, (14) 

respectively with d := (0,0,1). Equations ([^-(|^ and (j^ yield the following compatibility 
equations: 


r =^d, 
R' =RK, 


(15) 

(16) 


where K is the unique skew symmetric matrix satisfying k = axial (K) and R is the matrix of R 
relative to the fixed basis { 61 , 62 , 63 }. Henceforth, all components written with respect to the 
fixed {6,} basis are denoted by an over-tilde, e.g. r' = fje,-, f = {fi,f 2 ^rf), while components 
expressed with respect to the convected jd,} basis are written in san-serif font as in (11 1 . 

We employ the stored energy density according to the Kirchhoff model, viz.. 


W (k) = i (£/i xl+Eh 4 + GJ^ 3 ) 


(17) 


where E denotes the Young’s modulus of elasticity, G denotes the shear modulus, f is the 
area moment of inertia, £=1,2, and J denotes the (weighted) polar area moment of inertia 

m. 

We now normalize the system’s variables to non-dimensional form for a strip of length 
27r as follows: 


s = 2% —. 
L 


d _ 
ds 

n,- = 


2% d 
T df’ 

^%^Eh ’ 





m; = 


m,L 
2nEIi ’ 


(18) 

(19) 
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where L is the original length of the strip. With the exceptions of placing an over-bar above 
each quantity, the differential equations in ([T^-([T6|) are unchanged by this normalization. The 
normalized constitutive relations based on the stored energy density in 0 are 




fhi = ki, 


( 20 ) 



m 2 = A k -2 , 


( 21 ) 



m 3 = 7^3 , 


( 22 ) 

where 

A _ 

and 7 = 

GJ 

¥h' 

(23) 


Assume that the rod’s cross section is rectangular with width of length w in the di direction 
and and height of length h in the d 2 direction. The parameter A is 



for rectangular cross sections where vi> and h are the width and height of the cross section 
respectively normalized by the rod length L. Note that A = 1 corresponds to a rod with equal 
bending stiffnesses, while A<<lorA>>l corresponds to a rod with one very compliant 
bending direction and one very stiff bending direction, e.g. a thin strip. 

Formulas from strength of materials (cf. mi) give 


7 = 


2 

T+v ■ 


(25) 


where v is Poisson’s ratio. For numerical calculations, V = 1/3 is used which corresponds to 
7 = 3 / 2 . For clarity we now remove the overbars, with the understanding that all quantities 
are henceforth normalized according to 


3. Elastic Rod Solution Method 


Following the approach in ii, mi and we search for closed loop solutions of ([T3|)-([T^ 
that posses diflip symmetry about, say, the 62 , axis. That is, we suppose that a rotation by 180 
degrees of the closed rod about the 62 axis leaves the configuration unchanged. Hence we 
solve ([T3 ]i-([T6]i for half of the rod with appropriate boundary conditions (detailed in section 
|3.2| i, and generate a full loop solution by symmetry (detailed in section 3.4 1 . 

The resulting two-point boundary value problem for the half rod is solved on [0, n] using 
numerical continuation via the software package AUTO ifTSl . From this half solution, we 
generate a solution for the full Mobius strip on [0,27r]. As in ||7], the continuation is started 
from the equilibrium configuration of a twisted rod with equal bending stiffnesses, and then 
the path of equilibria is followed as the constitutive parameter A increases with 7 fixed. 


3.1. Parameterization 


Following the treatment in 13, ^ in ( [f^ is parameterized via quaternions, thus avoiding the 
usual singularities associated with Euler angles. Accordingly (15 1 -( 16 1 are replaced by 


f' = R (q) d, 
q'=A(q) k, 


(26) 

(27) 
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respectively, with q := {qo,qi,q 2 ,q 3 ) subject to the normalization 


^ 0 +? 1+?2+?3 — 1 ■ ( 28 ) 

The quaternion parameterization of the rotation matrix, ^ (q), and the quaternion differential 
equation matrix A (q) are (cf. El). 


/q'o + q'l-1/2 qiqi-qoqj 

= qiqi + qoqi ql+ql-^/'^ 

\q\q3-qoq2 qiqs+qoqi 


qm + qoqi \ 
qiq^-qoqi 
?o + 93-1/2/ 


(-qt 

-q 2 

-q3\ 

qo 

-q2 

q2 

q-i 

qa 

-qi 

\~q2 

qi 

qo / 


(29) 


(30) 


In general, an accurate numerical solution of (26l-(27l (satisfying reasonable boundary 
conditions) need not satisfy (28i with accuracy. We follow the approach in El and replace 
I with the augmented equation containing a multiplier /r S K: 


q'=A(q) k + ^q. 


(31) 


Use of ( [3 T] i ensures that (28i will be satisfied identically along the entire length of the rod 
whenever (|28|l is merely enforced on the boundary points. In practice, it turns out that the 
multiplier p takes on numerical values close to zero (typically p — O (lO^^)), cf. El . 
Combining 0 -([l4|), p0|)-(|22|), (j2^ and ( ^ we arrive at the full governing system: 


where 


n' + k(m) X n = 0, 

(32) 

m^ + k(m) X m + dx n = 0, 

(33) 

f'-^(q) d = 0, 

(34) 

q'-A(q) k(m)-^q = 0. 

(35) 

r t \ ( '”2 '«3\ 

k(m):= 

(36) 


3.2. Boundary Conditions 


Equations (32i-(35i constitute a system of first order ODE’s in 14 unknowns (f,n,m,q,/r) 
and the material parameters X, 7. Thus the problem requires 14 boundary conditions. We 
consider a half-rod of length n (so that the Mobius strip has total length 2n). The position and 
orientation of the rod at i = 0 are fixed at a particular point in yielding 


r(0 )-ei =fi (0) =0, 

(37) 

r (0) . 62 = f2 (0) = 1, 

(38) 

r(0)-e3 = f3(0) =0, 

(39) 

q(0) = {qQ,qi,q 2 ,q 3 ) = (1,0,0,0) . 

(40) 
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-0.5 


Figure 2. The rod in the initial equilibrium configuration used to start the continuation 
calculations. It has a fixed end at x = 0 and a hinged end free to slide along the e 2 axis at 
s = 7C. The coloring is presented to show orientation of the material points of the rod as the 
cross-sections rotate clockwise about d 3 from i = 0 to x = ;r. The complete closed rod is 
formed via a reflection about the 62 axis. 


This gives seven boundary conditions. We note that (40i fulfills (28 1 at the “left” boundary 
point. 

For the boundary conditions at s = ;r, we place a perfect hinge parallel to 62 , along which 
the rod may freely slide and about which freely rotate. From this, the boundary conditions at 
s = n become 


r(7t:)-ei (tt) =0, 

(41) 

r (tt) • 63 = f 3 (tt) = 0 , 

(42) 

n(7r)-62 = ni (tt) = 0 , 

(43) 

m(;r) • 62 = mi (;r) = 0 . 

(44) 


In addition, at the end the rod will twist by a quarter turn relative to the orientation in (40 1 . 
Thus the directors at i = tt have the form 


di {n) = -62, 

d 2 (tt) = cos j3 Cl + sin j3 63 , 

ds (tt) = — sinj3 Cl + cos j3 63 , 


(45) 

(46) 

(47) 


where p is some unspecified angle. Comparing the rotation matrix specified by (|43|)-(|4^ with 
the parameterization of R (q), cf. (29 1 , we may choose the two non-redundant conditions 


‘7o(^)+d(^)-^ =0, 

qi (?t) ( 7 t) - qn {n) q\{n) = 0 . 

In addition, we impose the normalization 

q\ {T^)+q 2 {T^)+ql {T^)+ql (?r) = i, 


(48) 

(49) 

(50) 


fulfilling (28 I at the boundary point. This completes the required set of 14 boundary 
conditions. 
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3.3. Initial Equilibrium Configuration 

As in Q, the starting equilibrium configuration for our continuation scheme is a rod with 
X = 1 (i.e. circular cross-sections made of an isotropic material) and 7 = 1.5 deformed in a 
semi-circular configuration which is rotated clockwise (from the viewpoint of s = 0 ) by a total 
angle of ;r/2 at s = ;r, (Figure |^. The configuration for the half-rod defines the configuration 
of the full Mobius strip on s € [0,2;r] via reflection about the 62 axis. 


3.4. Full Rod Construction 


Once a numerical solution for @-(|4|), (l4g-(igi, (l48|-@ is obtained for 

s G [ 0 , 7 r], the rest of the closed-loop for s G [K,2n] is constructed via a rotation by 180 
degrees about the e 2 -axis. The following procedure is rigorously detailed in || 8 ]. Denote the 
calculated solution to (32 1 -( 35 1 for s € [0, n] by a superscript “c”, e.g. (s) for the calculated 
rod centerline position. The flip across the axis of symmetry is given by 


E = - (ei (g)ei) -I- (62 062) - ( 63063 ) . 
The position of the centerline for s G [0,27r] is given by 

r^' (s) s G [0, n] 


r(s) = 


Er‘'( 2 a:-s) sG[n,2n\ 


( 51 ) 


(52) 


The continuity of r() dX s = K follows from ( |4T] i, and (|5^. It follows similarly that 
r(0) = r(2;r). 

The extension of the rod’s orientation on [7i:,27i:] is defined in terms of the director fields 

d;(i); 


dl(.) = j 

|d^(.) 

s G [0, %] 

(53) 

[Ed^ {27t-s) 

s G [n,2n] ’ 

d2 (s) = j 

f d| (s) 

[-Ed^ (27t-s) 

s G [0, n] 
s G [it,2n] ’ 

(54) 

d3 (i) = j 

idUs) 

[-Ed§ (27t-s) 

i G [0, n] 
s G [%,2%\ 

(55) 


Observe that d,- (•), i = 1,2,3 is continuous on [0,2;r]. 

While the rod position and orientation are enough to reproduce the equilibrium 
configuration for the closed loop, we also give the extensions of the contact force and contact 
couple fields, n and m respectively, mainly for use in the stability analysis presented in Section 
1^ Following the results in ISI, the required extensions are given by: 


n(s) = 


m(i) = 


n" {s) 

^-En''(27r-i) 
01 “^ (i) 

-Em'^ {2n-s) 


s G [0, 7Z] 
s G [7i:,27r] ’ 

i G [0, 7t] 
s G [%,2%\ 


(56) 


(57) 


In view of (43 i-(44 1 , we see that n (•) and m (•) are each continuous ats = n. We further claim 
that n (0) = n {2n) and m (0) = m (2;r). To see this, note that the global balance of forces and 
moments for the half rod on [0, ;r], together with (43 i-( 44 1 reveal that n (0) • 62 = m (0) • 62 = 0. 
The claim now follows directly from 
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(a) A = 1 


(b) A = 2 



(c) A = 10 


(d) A = 1000 





Figure 3. Kirchhoff rod theory results. All plots have y = 1.5. Note that the origin of 
the axes is at the point (—1.3, —1.3,0). Coloring is used to indicate the orientation 
of the strip. 


4. Kirchhoff Rod Theory Results 


The system (32i-([3^, subject to boundary conditions (|j7]i-(|40li, (41 i-(44i, and ( |48]l-(|50ll i s 
solved via collocation methods by AUTO-O7p Continuation and Bifurcation software ifTMlSl . 
The rod is divided into 30 mesh intervals with 4 collocation points per interval for a total of 
121 nodes along the interval [ 0 , 7 r], and the mesh is updated every three continuation steps. 
These numbers were chosen based on recommended values given in ITh), and increasing the 
number of nodes further did not lead to a significant quantitative difference in the numerical 
results. 

Starting from the equilibrium configuration with A = 1 and 7 = 1.5, new equilibrium 
configurations are found as A is increased and 7 is held fixed, cf. (36i. The calculated 


equilibrium configurations in Figure]^ are in qualitative agreement with the smoothly varying 
configurations found in ||2|. We note that different values of 7 in the allowed range do not 
produce qualitatively different equilibrium configurations. These updated results confirm that 
Kirchhoff theory does not capture the sharp localized bending and twisting seen in naia. 
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5. Local Stability of Kirchhoff-Model Configurations 


The local stability of the equilibrium configurations of the entire closed Mobius strip on [0,27r] 
is investigated in this section, for small but arbitrary perturbations - in particular, perturbations 
that break the flip symmetry of the equilibrium configuration. For any computed solution of 
the half-rod, the first step is to generate the entire solution on [0,27r] via (511-(57 1 . We then 
employ the methodology of ifTOl . The latter is quite general and can accommodate an accurate, 
discrete-numerical representation of a rod equilibrium in the presence of constraints - 
regardless of the numerical discretization method used, e.g., finite differences, finite elements, 
collocation, shooting methods, etc. Here we have the point-wise constraints inherent in (|^. 
We now summarize the methodology. 


5.7. Formulation 

We first introduce the spatially weak forms associated with the dynamics of unshearable, 
inextensible rods of length In without the presence of body forces or body couples: 


In 


Gd 


'ynamic 


pAr ■ p -l-p(Iw) • t/ 


ds, 


In 

Gstatic = j [n-(p'-v/xr')-l-m-v/' + ... 
0 


...-f <*„r'-Rea-f 1*3 (r'-Rea-l)] ds-f [n -p -fm-y/'lg'', 


(58) 


(59) 


where p is the density of the rod, A is the cross-sectional area, I is the moment of area tensor, w 
is the angular velocity of the cross-sections, the over dot (') indicates a derivative with respect 
to time, and p,i/, and ^ correspond to smooth variations in r,R, and n respectively. The 
spatially weak form of the governing partial differential equations governing the dynamics of 
the rod is represented by 


G (r, R, n , p, ^ ) — Gdynamic Gstatic — 0 5 (60) 

is satisfied identically at an equilibrium (r,R,n), for all smooth variations {p,\j/,^). We now 
consider small perturbations from an equilibrium (r, R, n) via 

re=r-)-eAr, (61) 

Re = exp (g A0) R, (62) 

ne=n-feAn, (63) 

where Ar, An are smooth admissible variations, exp( ) denotes the matrix exponential, A0 
is a smooth admissible skew-symmetric matrix, and e is a small parameter. We define 
A0 = axial (A0) along with the vector 

ACo=([Ar] [A0] [An])^. (64) 

The time dependent perturbations take the form 


AC=ACo^"'. 


( 65 ) 
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Taylor’s expansion about an equilibrium point generates 


G(re,R£,n,) = eDG(r,R,n)AC + o(e|AC|). 


( 66 ) 


Substituting (651 into the linear part of (|66|l, we obtain the generalized eigenvalue problem 


DGstatic ^Co — B ^^dynamic ^Co • 


(67) 


where p := is the eigenvalue. As discussed in ifTOll . the structure of (67i is nonstandard, 


due to the presence of the linearized constraints, e.g. on the left side, which are equated 
to zero on the right side. Moreover, for conservative problems, like the one at hand, the 
eigenvalues are necessarily real; A negative eigenvalue, p < 0, indicates instability, since 


(J = \/—p in (651 engenders exponential growth; a positive eigenvalue implies that CJ is 
purely imaginary, showing that ( [65] l is oscillatory. Accordingly, the solution is stable if all 
eigenvalues are positive. 

Explicit forms of DGstatic and DGdynamic for an unshearable-inextensible rod are derived 

inQol. 


5.2. Numerical Implementation 


To calculate the eigenvalues in (67 1 , we employ the finite-element method, as implemented 


in nniini. We approximate the smooth test functions (p, y/, ^) and spatial perturbations 
(Ar,A0,An) with piecewise linear functions. Nodal values for the variables (r,R,n,m) on 
[0,2;r] are obtained from the continuation results on [0, ;r] and symmetry transformations 
in (52i-(56i on [n^2n\. For N elements, this discretization transforms (67 1 into the matrix 
eigenvalue problem 


mxm x p 

^ 0 

pxm ^pxp 




Ar 






Ar 




A0 



X m ^ 

0 0 



A0 




Ano 






Ano 



( 68 ) 


where K is the global stiffness matrix, C is the global constraint matrix, M is the global mass 
matrix, m = 6N, and p = 3N. Note that p represents the total number of point-wise constraints 
on the discretized rod. Also note that K is block tri-diagonal and symmetric for conservative 
loadings at equilibrium ini. 


As mentioned before, the constraint terms cause (67 1 and ( 681 to be singular. The Q-R 
factorization of of C has the form 


C = [ Q1 Q2 


R1 

0 


= Q1R1, 


(69) 


Following the procedure in ifTOl . we generate a non-singular reduced version of ( |68l l via 

(Q2^ KQ2) ACo = P (Q2^ MQ2) ACo, (70) 

KACo = pMACo, (71) 


where K and M are the projected stiffness and mass matrices. This eliminates all the spurious 
eigenvalues and reduces the total dimension of the problem from m + pto m — p. 

The Mobius strip problem considered here is conservative, and the projected mass matrix 
M is positive definite. Accordingly, the latter may be replaced by the identity matrix without 
impacting the signs of the eigenvalues in (j6^, and our final form of the eigenvalue problem is 

KACo^pIACo, 


( 72 ) 
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Table 1. The four smallest eigenvalues of DGsmtic for a Kirchhoff rod with N = 240 elements 



Smallest 

2nd smallest 

3rd smallest 

4th smallest 

A = 10 

1.07e-4 

0.0449 

0.06895 

0.141 

A = 100 

0.00126 

0.0490 

0.0730 

0.163 

A = 1000 

0.0118 

0.0510 

0.0794 

0.169 


where positive eigenvalues of K indicate stability and negative eigenvalues indicate unstable 
perturbations. Since the problem at hand is conservative, K is the discreteHessian 
corresponding to the constrained potential energy. Thus, (0, while derived as part of a 
linearized stability method, the positivity of K is equivalent to the minimum-potential-energy 
criterion. 

5.3. Boundary Conditions 

For the closed loop, both the position and the orientation of the rod at s = 0 and s = 2n 
are clamped. Assuming the rod is divided into N elements with N +l nodes, the boundary 
conditions are 


Ar(‘’)=0, A0(‘’)=O, (73) 

Ar(A'+i)=o, A0(^+^)=O. (74) 

These ensure that the s = 0 and s = 2n ends of the rod will remain smoothly connected under 
any perturbation. In addition, eliminate the six neutrally stable rigid-body modes 

corresponding to uniform translation and rotation of the closed rod and also one additional 
neutral degeneracy associated with the axial motion of the strip acting through its own fixed 
configuration ifTSl ISll. For the initial isotropic configuration (A = 1), there is one remaining 
zero eigenvalue, due to a one-parameter family of equilibria corresponding to continuous 
precession of the centerline configuration accompanied by rolling of the cross sections in the 
opposite sense, cf. ||8l [18] [191 • This degeneracy disappears for A > 1. 

5.4. Results 

The numerical equilibrium solutions from AUTO calculated in Section [^ are extended to the 
full Mobius strip on [0,2;r] and used for the finite element calculation. This results in a mesh 
resolution of 240 elements for the full rod. We find the eigenvalues of K using the eigs() 
function in Matlab for each equilibrium configuration. 

We list the three smallest eigenvalues of the Hessian K in Table for several values of 
the the bending stiffness ratio “A”. In all cases the computed eigenvalues are positive, and we 
conclude that the Kirchhoff-rod equilibria for the Mobius strip are stable with respect to all 
local perturbations - symmetric and anti-symmetric. 

6. Developable Rod Model 

In this section, we model the Mobius strip as a developable surface as in min [5]. We employ 
the developable-strip plate model of Wunderlich |[3l as derived in IQ. Once equilibria are 
calculated, their stability can be assessed by adapting the approach of cni as employed in 
section 0 
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Using the same notation introduced in section]^ the developable strip is defined with 
centerline coordinate s G [0,L] in a straight, stress free reference conhguration. The position 
of the strip is dehned by the vector-valued function r (s) with the reference conhguration’s 
centerline given by ro (s) = sea. The directors are again dehned in Q-([^, and the dehnitions 
@-0 remain valid. 

Departing from Cosserat rod theory, we assume that the strip has an instantaneous axis 
of bending given by the vector, b(s) S spanjdijda}. Let (j) denote the angle between the 
centerline tangent vector, da, and the instantaneous axis of bending. Dehne the quantity 
j] = cot(j). Note that rj =0 corresponds to the usual Cosserat rod theory. 

Inherent in the approach of ID, echoed in Q and Q, is the tacit assumption that the hat, 
stress-free reference conhguration admits the representation 


X = sea-hv [ei-l-T] (s) ea] , (75) 

where s G [0,L], v G [—w/2,w/2]. That is, the mapping (v,s) —>■ (v,s-|-v77 (s)) should be 
locally injective on Q := [0,L] x [—w/2,w/2], viz., 1 -|-vT7'(s) > 0 on Assuming this is 
the case, then the deformation of the strip, f; > E^, is given by 

x = f(X) =r(s)-fvb(s), (76) 


with 


b(s) =di(s)-hT7(s) da (s) , 


(77) 


Note that (76 1 dehnes a ruled surface with normal vector N = da (s). 

The strip is presumed inextensible and the centerline in ( [76] ) 
inextensible and unshearable via (0. In addition, the constraint 


is constrained to be 


k-2 = 0, (78) 

precludes bending along the stiff axis. The ruled surface in ( f76| l, ( fTTj l is developable if (cf. 

El) 


ffa - T] fCi = 0. (79) 

The strip ( |7^ is a tangent developable with one generator of curvature, given by b(s). The 
constraints (0, ( |78j l, and ( [79| l enforce developability. 

Following 1^ the stored energy for the thin strip is derived from a constrained St Venant- 
Kirchhoff plate. In particular, the only contribution to the stored energy is that due to pure 
bending about the instantaneous axis (0; an integration across the width yields the total 
stored energy expression due to Wunderlich: 


L 



0 


where 


Eh^ 

12(1-v2) ’ 




^ f l+ri’w/2 \ 
rj'w \ 1 — rj'wjlJ ’ 


D:= 


( 81 ) 
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and where E is Young’s modulus, v is Poisson’s ratio, and Ki is defined in Q. Following 161 , 
( [80| l is now amended by integral terms involving the constraints Q, ( [78] l, and ( [79] l and the 
appropriate Lagrange multipliers: 

L L L 

U = V + J m 3 (k '3 — T] tfi) di + J m 2 k '2 di + J n • (r' — d 3 ) di (82) 

0 00 


In order to derive the Euler-Lagrange equilibrium equations, we consider smooth, L- 
periodic variations f, ©, fj, where 0 is skew-symmetric valued, as follows: 

r^-r-l-af, R ^ R-fexp (a0) R, Tj->Tj-|-afj (83) 

where a is a small parameter. We then find 

r'r'-fa (f'-fr'X 0) -l-o(a) , (84) 

K ^ K + aO'+ o{a) , as a->0, (85) 


where 9 := axial (0). We substitute (831-(85 1 into (82 1 , take the derivative of the resulting 
expression with respect to a, and then evaluate it at a = 0. A formal integration by parts then 
delivers the first variation condition: 


5U = 


f J ^ 

ri 1 r (irn') 

—, 

1 



+ 2DwKiri [l+J?^]g(wT7') + ... 


^ E. E 

... — m 3 K'i > Tjdi —y^n'-fdi — J {m'+ {d-i xn)} ■ dds = 0, (86) 


for all smooth variations f, 0 , fj, where the overdot, () indicates a derivative of a function with 
respect to the whole argument and prime, ()' denotes a derivative with respect to centerline 
arclength, s, and 


m 


1 =di •m:=DwK-i [l-fT]^] g (wr;') -1] m3, 


(87) 


Choosing fj =0 and 0 = 0 for all variations f in ( 861 yields (j^, and the fj =0 for all variations 
0 delivers (j^. A new, third Euler-Lagrange equations, corresponding to all variations fj in 
,reads 


d 

ds 


+ g{wri') 


+ 2DwKiri [l+ri^]g(wri') - m^Ki =0. ( 88 ) 


We observe that ( 881 is singular, corresponding to the pointwise vanishing of Ki. The latter. 


in view of (78 1 , is the total curvature of the center-line curve, i —r (i) and our director basis 
here coincides with the usual Erenet-Serret frame. As such, the validity of (|g, and ([8§ 
requires Ki f 0 ,cf. lISlfTTI. 

Our intended strategy here is the same used before in section]^ viz., solve the governing 
equations on [0,;r] and then generate the rest by rotation (flip symmetry). In particular, 
the latter, purely kinematical requirement (flip symmetry of the configuration) implies that 
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Ki (n) = 0. In order to overcome the otherwise certain numerical difficulties associated with 
that, we introduce the following “elliptic regularization” into the energy functional: 


U = U + 




(89) 


where e > 0 is a very small parameter and U is the new regularized energy. This, in turn, 
modifies the principal part of as follows 


■ + I I7"+T>w^k-i [1 + i?^]^^(wT7') k[+... (90) 

...+2Dw^Ki [l + 1 ]^] g(wT 7 ') rj rj' -2DwKf [l Tt]^] g{wr]') J) + k-i m 3 = 0 
We normalize all variables and the strip length to 27t: according to ([T8]l and 


2’Kw 


n, 


2% m; 


w=—— = 

L Dw DwL 

Note that J] is already a unitless parameter on [—1,1]). The complete system of differential 

equations for the developable rod (dropping all overbars) is given by 

n' + k(m) X n = 0, (92) 

m'+ k(m) X m + d X n = 0, (93) 

r'-R(q)d=0, (94) 

q'-A(q) k(m)-^q = 0, (95) 

02 [j?,i?'] i?" + oi T?'+ao [i?,!?'] =0, (96) 

where 

d:=( 0 , 0 ,l) 

' (mi+qma)^^ t] (mi+7] m3)L^ \ 


k(m):= 


^47t^ [I+ ri^] g(wri') 


0 , 


02 [17,17'] =e- 
01 [ 17 , 17 '] 




mi + 1 ] m 3 

fl+T72]g(wT7')' 


4;r2[l + Tj2j g(wri') 

^ f g(wT70g (w?70-2 g'(w?70 
g(w77') 


(97) 

(98) 


/^2' 


(m 3 [ 1 - 17 ^] - 2 miT 7 ) (mi+Tjm 3 )g'(wT 7 ') 


« n'] - (’T^I+^ ' 213 ) 

oo [t7,^ J ^ 2 


[I + T72] g{wri') 

g(w 77 ') m 3 [ 1 - 17 ^] - 2 T 7 mi 


[I + T 72 ] g(wri') 


W 02 - '-fr + 

g(wri') 


[1 + T 72 


(99) 

( 100 ) 

( 101 ) 


and where the overdot, (), indicates a derivative of a function with respect to the whole 
argument and the prime, () , denotes a derivative with respect to centerline arclength, s. 

From (81 1 , we see that the stored energy in (80 1 has a diverging integrand when rj'w = 0. 
This is a removable singularity which will occur at the boundary point s = 0 due to (121 1 . To 
avoid division by zero, g is expanded in the following Taylor Series about wrj' = 0: 


g(W) = 1 + 


(wTj')^ (wri'f (wri')° (wT]') 




r\6 




12 


80 


448 


2304 


+ 0 


([wT]'] 


10 


( 102 ) 


This Taylor series is used instead of the exact expression in (81 1 in the neighborhood of points 
where rj' =0. 
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7. Developable Rod Solution Method 


The system of ordinary differential equations in ([92|l-(96 1 is solved using the same procedure 
outlined in section]^ AUTO is again used to produce a solution to the half problem on [0, k] 
and the full Mobius strip on [0,2?!:] is constructed through a flip about an axis of symmetry. 
However, a major point of departure from section is that there is no explicit (twisted) 
solution from which to initiate continuation. Instead an appropriate initial configuration is 
found utilizing a series of intermediate continuation calculations. 


7.1. Initial Equilibrium Configuration 

To form a twisted Mobius strip, the centerline must undergo a non-zero twist, Kj, which 
induces non-zero curvature, Ki, and a non-zero scalar parameter, t], via the constraint ( |79] l. 
Thus, the most direct closed form equilibrium configuration is an untwisted strip with rj = 0 
which satisfies ( [90| ) trivially. Only a strip in the shape of a cylinder satisfies rj =0 while being 
developable lET so the initial configuration is a strip with a semi-circular centerline on [ 0 , 7t] 
with uniform curvature Ki =n. 

In order to sustain the uniform curvature in this configuration, a moment on the s = n 
boundary will be applied. This will initially violate the zero-moment condition for the hinge at 
i = TT, so we use continuation to follow the path of equilibria to a non-uniformly curved strip 
with zero moment at the s — n boundary (detailed in Figure a)-(b)). Once this is completed, 
the hinge boundary condition is fixed by setting mi (tt) = 0. Next, another continuation step 
is used to twist the strip by n/2 and align the hinge with the axis of symmetry, say 62 (detailed 
inFiguregc)-(f)). 

After both of these continuation steps are completed, a Mobius strip is formed via a 180 
degree rotation about the axis of symmetry. Continuation is performed in e to reduce the 
regularizing term followed by continuation in the width w to find the Mobius-strip equilibria. 
In summary: 

(i) Start with an untwisted strip with constant curvature. This requires an applied moment 
at s = ;r to sustain the curvature. 

(ii) Use continuation to relax the moment until there is no applied moment at s = ;r. 

(iii) Fix the moment at zero, and use continuation to twist the orientation at i = tt until the 
strip will form a Mobius half-strip. 

(iv) Perform continuation in the width and parameter e. 


7.2. Boundary Conditions 

As in Section [3^ the half-strip satisfies p7|)-(|44|) and (48l-(50i. In addition, we require 

i?'(0)=0, (103) 

r]{n) = 0, (104) 


which together with the other boundary conditions ensure flip symmetry of the complete 
Mobius band with the generator of curvature at s = n, b(;r), is aligned with the axis of 
symmetry 62 . 

As explained in Section [TTT] we require some preliminary continuations steps in order 
to arrive at a starting Mobius strip configuration. Each stage of the preliminary continuation 
calculations requires its own set of 16 boundary conditions as detailed below. 
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7.2.1. Initial Configuration When starting from the uniformly curved conhguration (see 
Figure]^ a)), the boundary conditions for the hxed end are dehned as above in ([J7]i-(|44|. As 
in the Kirchhoff rod case in section]^ only the displacements in the ei and 63 directions are 
constrained at s = ;r 


r(;r).ei=ri(;r)=0, (105) 

r( 7 r)-63 = f3 (tt) = 0, (106) 

Ats = 7Z the positions are held fixed via 

r( 7 r).ei =fi( 7 r)= 0 , (107) 

r(7r)-e2 = f2(7r) = -l, (108) 

r(;r )-63 = r 3 (;r) = 0, (109) 

while the hinge with variable end moment is given by 

ni(;r)=0, (110) 

mi(7r) = 7r(l-Si) , (111) 


where the continuation parameter Ei S [0,1] begins at 0 and is continued to 1. For the strip 
orientation, the director di is constrained to point along the hinge axis via 


qi{n)q2{n)+qo{n)q3{n)=0, (112) 

and the the quaternions are normalization through 

(^)+?2 (^)+?3 (^)+ 9'0 = 1 ■ (113) 

The hnal two boundary conditions are given by (|103|l-([T04|). 


7.2.2. Moment Relaxing Continuation At s = ;r the positions are held hxed via 


r(7r)-ei = n (tt) =0, (114) 

r(7r)-e2 = f2(7r) = -l, (115) 

r(7r)-e3 = f3(7r) =0, (116) 

while the hinge with variable end moment is given by 

ni(;r)=0, (117) 

mi(;r) = ;r(l-Si), (118) 


where the continuation parameter Ei € [0,1] begins at 0 and is continued to 1. For the strip 
orientation, the director di is constrained to point along the hinge axis via 


qi{n)q2[n)+qQ[K)q3{n)=0, (119) 

and the the quaternions are normalized through 

q\{K)+ql{n)+ql{K)+ql{n) = l. (120) 

The hnal two boundary conditions are 

rj'(0) = 0, t7(;r)=0, (121) 


which ensure the s = % end of the strip is aligned with the axis of symmetry 62 - 
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7.2.3. Twisting Continuation For the twisting continuation, the i = 0 end remains via ( |J7| )- 
(ig. As in the Kirchhoff rod case in section]^ only the displacements in the ei and 63 
directions are constrained at s = ;r 


r(7r) -ei = n (tt) = 0, 

( 122 ) 

r(7r)-e3 = r3(7r) = 0 , 

(123) 

while the hinge conditions are enforced via 


ni( 7 r)= 0 , mi( 7 r)= 0 , 

(124) 


where (124 1 is equivalent to fixing Ei = 1. The strip needs to be twisted by ;r/2 and the end 
oriented along the 63 axis of symmetry. This requires the transformations 


di(;r)^-e 2 , da (;r)- 62 ^ 0 . (125) 

The rotation that executes ( |125| l is facilitated by the dummy continuation parameter Z 2 , which 
starts at 0 and is continued to 1. Define 


^23 = [d3(7r)-e2](Sj=i S2=0) > (126) 

which corresponds to the 62 component of the da vector when the hinge is fully relaxed and 
the strip is still untwisted. The boundary conditions 


2 {q\ ( 7 t) q2{Tl)+ qt) (;r) ^3 (;r)) + S2 = 0, 

(127) 

2 {q 2 (n) qj, {tz) - qi {tz) qo (tz)) +^23 ^2 = 0, 

(128) 

q\ i^)+ql +ql = i > 

(129) 


execute the transformation (125i and twist the strip about its centerline by njl when E 2 
is continued from zero to one. The final two boundary conditions are again p21[ ). Once 
Z\—Z 2 = \, the strip is a half-Mobius strip on i S [0, n] and the boundary conditions are held 
fixed for continuation in the width, w and regularizing term e. 


7.3. Full Strip Construction 


Following the procedure outlined in section 3.4 and 13, the numerical solution is obtained for 
s G [ 0 , n], and the rest of the closed strip is constructed via a flip rotation by 180 degrees about 
the e 2 -axis. The transformations (51 i-(|57]i are applied along with the transformation 


Tj(s) = 


T7"(s) 

-ri^{2n-s) 


s e [ 0 , n] 
s e [n,2n] ’ 


(130) 


where the superscript c indicates the calculated value on the domain [0,;r]. Note that at 
s = 2n, the strip is twisted by K about the centerline, so dj {2n) is in the opposite direction of 
di (0). The negative sign in (1301 ensures that the strip forms a smooth closed loop after this 
orientation change. 


8. Developable Rod Results 

Due to the presence of constraints ( [78| ) and a finer mesh is needed to converge to 
equilibrium configurations for the developable rod than the one used in Section for the 
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©3 


©3 


©3 




(d) Twist continuation, (j> = 571/6 


(e) Twist continuation, (p = 27cl3 


(f) End of continuation, (p = 0, 
fonning half of a Mobius strip 


Figure 4. Initial and intermediate configurations to set up the Mobius strip. 0 is the angle 
between the d| (tt) and the axis of symmetry 62 - (a)-(b) Continuation in Sj to relax the mi 
moment about the hinge to zero, (c)-(f) Twisting continuation in S 2 to align the strip end with 
the axis of symmetry. Black lines indicate the direction of the generator of curvature, b (i). 


Kirchhoff rod. In particular, this ensures convergence during the twisting continuations steps 
detailed in Section 7.2.3 when 77 , 77 ' 0. The rod is divided into 100 mesh intervals with 5 

collocation points on each element for a total of 501 points for the rod on y G [0, tt]. The same 
continuation step size, tolerances and plotting routine from section|^are used. 

Starting with a hxed width, w, and regularizing term, e, continuation in the auxiliary 
variables, Ei and Z 2 , is used to arrive at the Mobius strip conhguration, as depicted in Figure 
1^ Once the Mobius conhguration is reached, continuation in the width w and the regularizing 
coefficient e are used to obtain the conhgurations shown in Figure]^ In particular, for a given 
hxed width, e is reduced as small as possible. However, it cannot be decreased indehnitely: 
As shown in Figures | 6 | 8 ] as e is decreased for given hxed width, the product W 77 ' approaches 
the value —2 at the hinge location s = n at which the density function ( [ST) blows up. As 
discussed after equation |7^ , |w 77 '| = 2 indicates a breakdown of injectivity for the mapping 
( |75l l at the edge of the strip, and the through-width integration leading to ( (80| l, ( [ST] ) in lO is 
no longer valid. Of course the right side of (81 1 is undehned for |w77'| > 2 We further observe 
from Figure that the larger the width, the larger the value of regularizing coefficient. In 
Figures and , we plot the computed “generators of curvature” on the reference strip, viz., the 
right side of ( |75| l at regularly spaced values of s with v ranging over the width. Observe that 
as e is decreased or w increased, the generators almost intersect at the midpoint (s = n) where 
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the hinge boundary condition is enforced. 

In Figurej^we illustrate some key solution helds for various hxed widths (as indicated). 

In each held we also demonstrate the robustness of our results for three consecutively 
decreasing values of the regularizing parameter. For each of the widths w = 0.2, 0.4, and 1.0, 
the plots for Ki and rj demonstrate singular-perturbation behavior in a small neighborhood 
of the hinge location s = 7Z. For the larger widths w = 1.6 and w = 2.0, this effect is less 
concentrated, which is undoubtedly due to the large values of e required. 


9. Developable Rod Stability 

As in the case of the Kirchhoff rod in section the local stability of the conhgurations 
from section is determined through linearized dynamic analysis about the equilibrium 
conhgurations via an adaptation of the procedure in IfTOl . 


9.1. Derivation of Stability Equations 


As in section HD the spatial weak form of the dynamical equations governing the rod-strip 
|-( 101 1 is given by 


2k 


G = 


(fr 


pA-^-n ') •p-l-(H-m'-r'xn) -1/4-... 


...-f ^3 (r'-da- 1 ) -l-atiK-- (da - qdi) +(02K-d2 + ... 


...+ 


-erj - 




2n^w d 
Ef ds - 

- + ^ (1 + 8 {wr]') - ma fCi 


(131) 


X 


where p is the density of the rod, A is the cross-sectional area, H is the cross-sectional angular 
momentum, (Oi, CO 2 , and x correspond to smooth variations in r,R, n, m 2 , ma, 

and T] respectively. As in the Kirchhoff rod case, H and the dynamic pieces are not needed 
to calculate the stability of a conservative system, so they are not explicitly derived. After 
integrating by parts, G decomposes into dynamic, static, and boundary contributions with the 
static piece given by 


2n 


Gstatic= / < n-(p'-v/xr') -pm- v/'-t-<^/r'-d/-f ^3 (r'-da-l) -f tUiK--(da - pdi)-f... 


, , , 2n^w 

... + ® 2 K-• da-f e p'x'-f 




Z' 


(132) 


...-f 


^ K-^T] (1 -f g (w 77 ') - ma K-i 


z 


The contact couple, mi, is constitutively determined, and hence, static solutions 
corresponding to the zeros of (|132|i satisfy the dynamical, spatial weak form 


G(r,R,n,m2,ma,77 ; p,v/,^,®2,®3,Z) =0, 


(133) 
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(a)w = 0.2 e = 0.033 


(b)w = 0.4 e = 0.11 




(c) w = 1 e = 1.0 


(d)w=1.6 e=4.08 



(e) w = 2.0 e = 10 

Figure 5. Mobius strip configurations for different widths. For each w value, the configuration 
with the smallest e value is shown. Note that the width in 0 corresponds to one-half the width 
in this formulation. 


Computation of Unconstrained Elastic Equilibria of Complete Mobius Bands and their Stability23 



. w = 0.4 e = 0.25 

- w = 0.4 0.1667 

-w = 0.4 e = 0.1111 





. w = 0.4 e = 0.25 

- w = 0.4 e = 0.1667 

-w = 0.4 e = 0.1111 



= 0.2 e=0.25 
= 0.2 e=0.1 
= 0.2 £= 0.0333 



. w = 0.4 £= 0.25 

- w = 0.4 e=0.1667 

-w = 0.4 e = 0.1111 


ature, twist, moment, rj, and T]fw values for different continuations 
configurations for w — 0.2 and w — 0.4. Note that rj shows a boundary layer around the 


point s — 71. Also note that the energy density in (ID becomes complex for w — ±2. 
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. w= 1.6 e= 5. 

- w= 1.6 e = 4.5 

-w= 1.6 e = 4.0816 





. w= 1.6 e = 5. 

-w= 1.6 e = 4.5 

- w= 1.6 e = 4.0816 




Figure 7. The curvature, twist, moment, rj, and rjfw values for different continuations 
configurations for w = 1.0 and w — 1.6. Note that rj shows a boundary layer around the 


point s — 71. Also note that the energy density in (ID becomes complex for w — ±2. 
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. w = 2. £=^2. 

- w = 2. e= 11. 

- w = 2. e= 10. 



. w = 2. e=^2. 

- w = 2. e= 11. 

-w = 2. e= 10. 


wq' 



Figure 8. The curvature, twist, moment, 7], and r]/w values for different continuations 
configurations for w — 2.0. Note that rj shows a boundary layer around the point s — 71. Also 
note that the energy density in (8T) becomes complex for w — ±2. 



(a) w = 0.4, £ = 0.25 



(b) w = 0.4, £ = 0.167 



(c) vy = 0.4,£ = 0.111 


Figure 9. Plots of the generator of curvature, b, on the Mobius strip reference configuration 
for w — 0.4 and different values of e, the regularizing term. The left end corresponds to 5 = 0 
and the right end corresponds to 5 = 2;r. 
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(a) w = 0.2, e = 0.033 



(b) w = 0.4,8 = 0.11 





Figure 10. Plots of the generator of curvature, b, on the Mobius strip reference configuration 
for different values of w. The left end corresponds to s = 0 and the right end corresponds to 
s = lit. 


at a static equilibrium. The linear perturbations 

f = r + aAr, R = exp(aA0)R, 17 = 77 +a At], 


(134) 


are applied where Ar, A0, A 77 are smooth admissible variations, exp( ) denotes the matrix 
exponential, and A0 is a smooth admissible skew-symmetric matrix. We define the vector 
A0 = axial (A0) and note that ( 134i induces the variations 

d; = d; + a (A0 X d,) , k,- = k-,- + a [A0' • di] . (135) 


In addition, the Lagrange multiplier fields {ni,n2,n3} and {m2,m3}, are linearly 
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perturbed by the expressions 


n = (n,'+ aAn,) (d; + a [A0 X d,]) , ; = 1,2,3, (136) 

rhf = (m^ +aAmf) (d^ + a [A0 X d^]) , i = 2,3. (137) 

We define the vector of unknowns, Ai^q and their time dependent perturbations, A^, as 

ACo=([Ar] [A0] At] [An] [Am])^ , AC=ACoe'"‘, (138) 

As in section [^Taylor’s expansion about an equilibrium point and substitution of (138 1 into 
the linear part produces the generalized eigenvalue problem 


C^Gstatic^^Q — Cr DG^yfiamic 0 • 


(139) 


where p := is the eigenvalue. As discussed in section]^ a negative eigenvalue, p <0, 
indicates instability, and the solution is stable if all eigenvalues are positive. 

Since the problem is conservative, the explicit form for DGdynamic is not needed, so it 
is not derived here. The derivation of the explicit form for DGstatic is provided in appendix 
[Appendix A The explicit form of DGstatic is 


Itc 


DGs, 


■ = j RAn ■ [p' — \ifx r') + (n x Ar') ■y/ — {nx A0) • (p' — x r')...+ 

v+... 


R 


^ 8 {wr\') R^ [A0'] • ei - t] Ama 
Ama 
Ama 


...+ 


4n:^wK'i (l+ 7]^)^ , 


L 2 


g' (wT]') Ap'R^’v/'-d + ... 


...+ 


IfiTT^fCit? (l + t?^)g(w77') \ / 

-1-ma ) ApRV- 


ei + A0 X R 


mi 

m 2 

ma 


(140) 

v+... 


...+R^-(Ar'+[r'xA0])+R 

167t:2 KiP (l + g (wT]') 


■ -0)1 7] ■ 


0)27] 

0)2 

•A0'+ k-iR 

0)1 (l-T]^) 



-0)2 


•A0 




L 2 


-ma ^ R^ [A0'] -ei-;i;ffiAma + e;i;'AT 7 ' + ... 


+ 8;r^wx ^2^ (l + Tj^)g(wT7') Atj' + ^^k-? {l + 3p^) g {wp') Ap + ... 


L 2 


47t:2x'wK'i fW , 2\2 ..t i\ s I 

L2- 12 (1 + ^ ) g{wp)Ap +... 

... + 2ki p {l+p^)g{wp') Ap + (1 + p^)^ g (wp') R'^ [A0'] -eilds. 
In addition, the boundary terms are given by 


- 0)1 KiAp 


Gbdry — 


n-p +myr + ep' z 


2n^w r 2 


L 2 


^•^(l + J?^) g{wp') 


In 


(141) 
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9.2. Numerical Implementation 


As with the Kirchhoff rod in section!^ the eigenvalues in (1391 are calculated via the finite 
element method as implemented in lUOl[TTll . The smooth test functions {p,\j/,^,co,x) and 
spatial perturbations (Ar, A0, An, Am, Atj) are approximated with piecewise linear functions. 
Nodal values for the variables (r,R,n,m,Tj) are obtained from the continuation results 
on [0,;r] and symmetry transformations detailed in section 9.3 For N elements, this 


discretization transforms (139 1 into the matrix eigenvalue problem 


Kmxm ^mxp 


CT 0 

'-'pxm ”pxp 



Ar 



A0 



1 

> > 

_ 1 _ 



Am 





M. 


mxm ^ 

0 0 



Ar 



A0 



> > 

_1_ 



Am 



(142) 


where K is the global stiffness matrix, C is the global constraint matrix, M is the global mass 
matrix, m = IN, and p = 5N. Note that p represents the total number of point-wise constraints 
acting on the discretized rod. 

This problem has the same structure as the stability calculation for the Kirchhoff rod in 
section]^ and it again results in the eigenvalue problem ( |72] l. Employing the same solution 
procedure from Qol , we find the eigenvalues of K where positive eigenvalues indicate stability 
and negative eigenvalues indicate unstable perturbations. 


9.3. Full Strip Construction and Boundary Conditions 


For the stability calculation, the closed developable strip is generated from AUTO’s half 
solution differently than the Kirchhoff rod case. If a solution on [0, n] was extended via 


(52i-(57 1 to a solution on [0,2;r], a node would be placed exactly at the s = 7t singular point. 
Instead, we extend the solution on [0,;r] to a solution on [—n,n\ via a flip rotation by 180 
degrees about the e 2 -axis to avoid placing a node exactly at this singular point. As in the 
Kirchhoff rod case, the following procedure is rigorously detailed in ISl . 

Denote the calculated solution to (921-(96 1 for s G [0, n] by a superscript “c”, e.g. (s) 
for the calculated rod centerline position. The position of the centerline for s G [0,2?!:] is given 
by 


r{s) = ■ 

/Er"(-s) 

|r^(5) 

s G [-^,0] 
s e [0, 7t] 

(143) 

where E is defined in ([sT]). The extension of the rod’s orientation on | 
terms of the director fields d; {s): 

— 7r,0] is defined in 

di (s) = 1 

\-Ed‘ii-s) 

s G [-7r,0] 
i G [0, 7t] 

(144) 

d2 (i) = 1 

fEd5(-s) 

[d^ (s) 

s G [-7t,0] 
s G [0, 7t] 

(145) 

da (s) = 1 

f-Ed^3(-s) 
[d^ (^) 

s G [-;r,0] 
i G [0, n] 

(146) 
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Observe that d, (-), i = 1,2,3 is continuous on [—n,7t]. Following the results in fSl, the 
required extensions are given by; 

(147) 

(148) 

In addition, the transformation for the scalar parameter t] and its derivative are 



s e [-?t,o] 
s e [0, 7t] 

s G [-7i:,0] 
i G [0, n] 


J](s) 


7]'(s) 


ri^i-s) iG[-7r,0] 

ri^is) sG[0,;r] 

-(jjr(-s) 5G[-;r,0] 


(149) 

(150) 


In our closed strip, both the position and the orientation of the rod at s = — ;r and s = 7t 
are clamped. Assuming the rod is divided into N elements with N +l nodes, the boundary 
conditions are 


Ar(“(=0, A0(“)=O, Ari^°'>=Q, (151) 

Ar(^+i)=o, A0(^+^)=O, At7(^+'(=0. (152) 

These 14 boundary conditions ensure that the s = —n and s = n ends of the rod will remain 
smoothly connected under any perturbation and satisfy the variation of ( |141| l. As before in 
Section [53| ( |151| l-( [T52| eliminate the six neutrally stable rigid-body modes corresponding to 
uniform translation and rotation of the closed rod and the additional degeneracy associated 
with axial motion of the strip acting through its own hxed conhguration. 

9.4. Results and Summary 

The numerical equilibrium solutions from AUTO calculated in Section]^ are extended to the 
full Mobius strip on [0,27r] and used for the hnite element calculation. This results in a mesh 
resolution of 1000 elements for the full strip. As shown in Table all the eigenvalues of 
DGstatic positive. Observe that the smallest (positive) eigenvalue consistently increases as 
the regularizing parameter is made as small as possible. We conclude that the developable-rod 
conhgurations are stable with respect to all sufficiently small perturbations -symmetric and 
non-symmetric. 

10. Concluding Remarks 

Here we present the hrst evidence for the local stability of flip-symmetric conhgurations of 
complete elastic Mobius bands. We employ two distinct models - the Kirchhoff rod model 
and the developable-surface model of Wunderlich. For the latter we present a novel strategy 
for the computation of complete unconstrained Mobius bands. To the best of our knowledge, 
it is the only known method for accomplishing such. Our introduction of a small elliptic 
regularization, cf. ( [89| , similar to what is often done in phase-transition problems, avoids the 
inevitable singularity associated with the rod-like formulation 13, lb) for Mobius bands, cf. 

Q, ca. 
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Table 2. The four smallest eigenvalues of DG static for a developable rod with N = 1000 
elements. 



Smallest 

2nd smallest 

3rd smallest 

4th smallest 

w — 0.2 

e = 0.25 

0.0073 

0.0335 

0.0409 

0.0589 

w = 0.2 

e = 0.1 

0.0079 

0.0359 

0.0376 

0.0508 

w = 0.2 

e = 0.033 

0.0087 

0.0329 

0.0339 

0.0483 

w — 0.4 

e = 0.25 

0.0077 

0.035 

0.0415 

0.0608 

w — 0.4 

e = 0.167 

0.0082 

0.0372 

0.0403 

0.0577 

w — 0.4 

e = 0.111 

0.0089 

0.0391 

0.0396 

0.0574 

w — 1 

e= 1.333 

0.0099 

0.0365 

0.0625 

0.109 

w = 1 

e= 1.125 

0.0102 

0.0386 

0.0603 

0.1049 

w — 1 

£ ^ 1 

0.0104 

0.0406 

0.0587 

0.1019 

w = 1.6 

£ ^ 5 

0.0164 

0.0412 

0.0928 

0.1681 

w = 1.6 

e = 4.5 

0.0167 

0.0436 

0.0907 

0.1643 

w = 1.6 

e = 4.082 

0.0173 

0.0478 

0.0875 

0.1573 

w — 2 

e = 12 

0.0245 

0.0484 

0.1254 

0.2206 

w — 2 

e = 11 

0.0245 

0.0508 

0.1216 

0.2147 

w — 2 

e = 10 

0.0255 

0.0575 

0.1155 

0.2008 


Unlike the numerical approach discussed in Q, ours presented here delivers complete- 
loop configurations in the absence of extraneous external fields. Moreover the solutions 
presented in Figures|6]|^and the eigenvalue results in Table|^all demonstrate the robustness of 
our results in the small regularizing parameter e. Finally we mention that our implementation 
of the Wunderlich model, for both computing equilibria and assessing their stability, is 
applicable to many other thin-strip problems. 
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Appendix A. Derivation of DGstatic for a Developable rod 

This appendix contains the work for the derivation of DGstatic in ( |140| i for a strip of length L. 
The quantity Gstatic is given by 


Gstatic — 


n- (p' - v/x r') + m- v/' + ^/r' -d; + 1*3 (r' - da - l) + ■■■ 


... + 0)1 K-- (da - T7di) + 0)2K--d2- 


... +eri X + 


w 


2 L 2 . 


+ g(wT7') 


Z' 


...+ 


2 ^ 

^Jfl2Tj(l+Tj2)g(wTj')-m3K-i 


Z M*. 


(A.l) 


The quantity DGstatic is the linearization of Gstatic about an equilibrium point, viz. 


DGstatic = ^ [Gstatic (t + ttAr,exp (aA0) R, n + aAn,m + aAm, rj + aArj ) ] (A.2) 

















Computation of Unconstrained Elastic Equilibria of Complete Mobius Bands and their Stability!) 1 


Using the notation 5{-} = [■]a= 0 ’ derive DGstatic term by term starting with the hrst 

component in ( |A.1[ ): 

5{n - (p'-v/xr')} = ^ [(Pi + aAn,) (d; + a [A0 x d;]) -{p'-y/x (r' + aAr)) ] 

(A.3) 

= An,d,' ■ (p'-y/x r') - n,d; • (y/ x Ar') + n; (A0 x d,) - (p'-y/x r') 

(A.4) 

= R An ■ (p' — y/x r') + (n x Ar') ■ y/— (nx A9) ■ (p' — y/xf) 

(A.5) 


The second term in ( |A.l| i is significantly more complicated because it involves quantities 
subject to variations in the rod orientation. The variation of the second term starts as 


= — [mid~i+m2d2 + m3d3]^^Q-v/' 


(A. 6 ) 



0 


0 

R 

A m 2 

+ A0 xR 

m 2 


A m 3 


m 3 


. Since m 2 and m 3 are both Lagrange multipliers, their variations are straightforward yielding 
^ [m2d2 + m3d3]„^o = ^[(m2 + aAm2) (d2 + a[A0xd2])]„=o + - 

^ [(m 3 + aAm 3 ) (d3+ a [A0 x d3])]„^o (A.7) 


(A. 8 ) 


For the variation in mi, the moment is constitutively determined, so it is significantly more 
complicated. From (36 1 m 1 is 

- 1 2 

fhidi = ^ (l + p^) g(wfi') kidi - r) m 3 di (A.9) 

= + (-(jj + a [A0 X di]) (^1 + [tj + g (w [tj' + aArj /]) + ... 

... - [tj + aArj] [m 3 + aAm 3 ] [di + a (A0 x di)] (A.IO) 
Note the vector quantities of the variations. Taking the derivative with respect to a yields 


da 


[™l3l]a=o = 


(l+TJ^) 


^2 g(wTj') [A0'-di] + 
4Tfi TJ (1 +TJ^) 


w 


K -1 (1 + Tj 2 
L 2 


-g'(wTj') Atj' + ... 


^2 S {wri') - m3 

J K -1 (1 TTJ^''^ 


L 2 


Atj — tj Am31 di +... 
g(wTj')-m3Tj i [A0 Xdi] . (A.ll) 


Since [A0'] is written with respect to the fixed basis we use 


[A0'-di] =R^ [A0'] -ei 


(A.12) 
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Adding this piece to ( |A.l 1[ ) and combining with ( |A.8| l, the total variation is 

5 {m- V/'} = 

^^^,2 ^ g {wrj') [A0'] • ei - tj Ama 


.+R 


.+R 


Am 2 

Ama 


^ g' (WT]') At]' + 


■w +■■■ 


- y -^g(w77')-ma 

0 

0 


At] 


...+A0 xR 


,2 gM ) 


^2 g(wT7')-maTj 

m 2 

m 3 


V + ... 


I/'. (A.13) 


Substituting ( [3^ into the cross product term and simplifying yields 
5{m- v/^} = 

^^^^2 ^ ^ (m'J?') R^ [A0'] • ei - tj Ama 

...+R ■r' + ... 

Ama 


.+R 




'^'^^^^(wTjO-ma 

0 

0 


Atj 


■ V^' + ... 


...+A0 xR 


mi 

m 2 

m 3 


■V- 


(A. 14) 


The third term of (|A.l|i is relatively straightforward and yields 


5{^,r'-d/ + ^a (*''■‘*3-1)} = ;^ 1 (r'+ aAr') • (di+a [A0 x di ])] + ... 


[^ 2 (r' + aAr')-(d 2 + a[A 0 xd 2 ])]„^o + ... 
•■•+d^ [^3 {(r' + «Ar') • (da + a [A0 X da]) - 1} ] . 


= ^,.Ar'.d,- + ^,.r'.[A0xd,-], 

= R^-(Ar'+[r'xA0]) . 

The second line, or fourth term, of ( |A.l| i is 

5[®iK--(da- T7di) + a)2K--d2] = 

- OTi (K' + aA0') • {rj + aArj) (di +a [A0 x di]) +... 

... + 0)1 [K + aA9') • (da + a [A0 X da]) + ... 

... + 0)2 (K' + aA0') • (d 2 + a [A0 x da]) . 


(A.15) 
(A.16) 
(A. 17) 


(A.18) 
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5[®iK--(d3- T7di) + a) 2 K--d 2 ] = 

-a)ifc-A77di +A0'- (-®i77di + ® 2 d 2 + ©ids) +... 

... + ®iff-([A0 xd 3]-77 [A0 xdi]) + ... 

... + ® 2 K--[A0 xd 2 ] (A.19) 


5[a)iK--(d3- Tjdi) + a) 2 K--d 2 ] = 

A0 • (-®2K'id3 + a)ifCid2)+A0 • (-a)i fCiT7^d2 + a)2i? k-i di) 
Putting into matrix form 

5[a)iK--(d3- 77di) + C 02 K--d 2 ] = 



■ -0)17] ■ 


®2l? 

0)1 ifi Atj +R 

0)2 

•A0'+ (fiR 

0)1 


0)1 


-0)2 


(A.20) 


(A.21) 


The variation of the Euler-Lagrange equation for the terms x’ is 


'{eri'x'- 


w 
2L^ L 


wx! d 

lEf da 




K^{l+ri^f g{wri’) x'] =£X'^ri' + - (A- 22 ) 

a [A0' • di] (^1 + a [t] + Aq]^^ g (w [t]' + aArj'] ) 


a=0 


= ex'ATj' + ^-^|(l+Tj2)^^(wT7')R^ [A0'] -611 + ... 

+ '^1 (l + J?^)^l(wJ?')A77' + 2K-iT7(l+T72)g(wT7') At 7 | (A.23) 

Finally, the variation of the Euler-Lagrange equation proportional to X is 


L 2 


Kjri (l+r^^) g{wri')-miKi 


X} = 


+ «Am 3 ) (tfi + a [A0' • di]) ] ^^^ +... 


da 

L 2 (Jo; 


(A.24) 


'^1 


+ a [A0'-di])^(T7 TaAt]) (^l + [rj + aArjf'j g (w [q' + aAq']) 


- a=0 


UkiT] {1 + ri^) g{wri') l„rrAo/i 
= Xi -^-m3 ^ [A 0 ] -ei -XK-iAm3 + ... 

+ (l+T7 ^)g(wT7 ') Atj' + ^K-? (l-f3TJ^)g(wT7') Atj (A.25) 

Combining the six terms in (|A.5| l, ( |A.14| i, (|A.17 | i,( |A.21| i, ( |A.23| i, and ( |A.25| l yields the 
formula for DGstatic given in (|140|l in section [9!l[ 
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